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ABSTRACT 

A simple and efficient method of artificial compression is introduced. This method is 
based on a modification of the slopes of the ENO reconstruction and, with the help of 
suitably chosen parameters, greatly improves the resolution of the contact discontinuities. 
Numerical examples are provided to test the performance of the method and to give some 
suggestions as to the choice of the parameters. 
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§1 Introduction 


In this paper, we introduce a simple yet effective method to improve the performance of 
the ENO schemes for hyperbolic conservation laws. 

Consider the following initial value problem: 

f 

( 1 . 1 ) < u < + /( u )* = 0 

u(x, 0) = u 0 (x) 

where 


u : R 2 R d 
/: R d — ► R rf 


Assume that the system (1.1) is strictly hyperbolic in the sense that the Jacobian 

( 1 - 2 ) a=^L 

du 

has only real eigenvalues with a complete set of eigenvectors. 

It is well known that the solution of (1.1) may develop discontinuities in finite time even 
though the initial value u 0 (x) is very smooth, say, a C°° function. These discontinuities 
include shocks, contact discontinuities and the wave fronts of rarefaction waves. 

It is natural that efforts in numerically simulating the solution of (1.1) with these struc- 
tures mainly focus on designing numerical schemes with the following properties: 
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i) achieving high accuracy in smooth regions of the solutions, 

ii) producing sharp profiles for the shocks and the contact discontinuities, 

iii) avoiding superfluous oscillations in front of or behind the shocks and the contact 
discontinuities, 

iv) getting correct positions and speed of the discontinuities, and 

v) avoiding non-physical discontinuities (e.g., expansion shocks). 

The Lax-Wendroff theorem (see [13]) says that, for a convergent conservative scheme 
consistent with (1.1), the limit function of the numerical solution as the mesh size tends to 
zero is a weak solution of (1.1). Thus, conservative schemes automatically guarantees iv). 
Throughout this paper, we shall only consider schemes in conservation form. 

To enforce v), one has to consider the so called entropy conditions which distinguish the 
physical solution from others. However, we are not going to discuss this problem here. 

Recently, Harten, Osher, Engquist and Chakravarthy (see [2], [3] and [4]) introduced a 
class of essentially non-oscillatory (ENO) schemes which are of globally high order of accuracy 
where the solution is smooth. Although a complete theoretical analysis of these schemes has 
not been done, extensive numerical examples show that these schemes are nonlinearly stable. 
Hence, these schemes excellently meet the requirements i) and iii). For ii), the ENO schemes 
produce extremely sharp shock profiles. However, they smear the contact discontinuities at 
a rate which appears to be of order 0(n r + 1 ), where r is the order of accuracy, and n is the 
number of the time steps. In order to overcome this difficulty, Harten introduced the concept 
of subcell resolution. This led to excellent results in 1-d computation. He is currently 
considering 2-d extensions. Another method which can be used to sharpen the contact 
discontinuity is Mao’s method introduced in [11]. Merging his ideas with subcell resolution, 
Mao’s method, hopefully, could be used with any known scheme. Some remarkable results 
have been obtained. However, much work needs to be done in order to make it practical in 
the cases of systems and of multi-dimensions. 

In this paper we will introduce a procedure which is different from Harten’s and Mao’s 
and which greatly improves the performance of the ENO schemes at the contact disconti- 
nuities. This is going to be achieved by combining the ENO schemes with a ACM (artificial 
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compression method, see Harten’s pioneering work [1]) type technique. In this way, we 
can effectively prevent the contact discontinuities from being smeared and, in most cases, 
limit the transition of a numerically computed contact discontinuity to about 2 cells while 
essentially keeping the order of accuracy in the smooth regions of the solutions. 

This procedure was originally invented for the cell average framework of the ENO schemes. 
Since then, Shu and Osher have transformed it to their pointwise ENO schemes. They have 
obtained surprisingly good results in 2-d computations. See [7] for the details. 

In the next section, we will give a brief review of the ENO schemes. §3 describes in detail 
our ACM type technique - the slope modification method and proves some of its properties. 
§4 applies the method to the system of Euler equations for gas dynamics. In §5, we will 
present some numerical results showing the performance of our method as well as giving 
some suggestions as to the choice of the parameters concerned. 

§2 Review of the ENO schemes 

We refer to [2], [3], [4], [5], [6] , [7] and [8] for details of the cell average ENO schemes 
and their pointwise versions. Here we only give a very brief review of the schemes. This is 
necessary for us to describe our slope modifying method. 

Originally, the ENO schemes were introduced for cell average values of the solutions as 
high order extensions of the Godunov scheme and the MUSCL scheme. Denote by y the 
numerical solution approximating the sliding average u of the exact solution u of (1.1), i.e., 

v(x,t) = u(x,f) 

1 r x + 

= hl^H U(V ' t)iy - 

Let X (<n,fn+i]} where x„ = ah, l„ = nr, j = 0,±1,±2,---, n = 0,1,2, be 

a partition of R x R + . Denote u(x it t n ) by uj and v(xj, t„) by Assume that E(t) is the 
exact solution operator of (1.1). I.e., if u(x,t ) is the solution of (1.1), then 

u(x, t ) = £(t)u(x,0). 
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The theoretical ENO scheme can be written as 

(2.1) v’ l+1 = i E{t)R(x, v n )dx 

h J X] _i 

which is said to be ‘theoretical’ since E(t) is exact. In (2.1), R(x, v n ) is a pointwise approxi- 
mation of the solution u at time t n derived from the cell average approximate solution v n . It 
is a piecewise polynomial. In fact, it is a polynomial in each cell (Xj_^, x J+ ^). The procedure 
of computing R(x , v n ) from v n is called reconstruction. 

Given the cell averages of a function u(x), in [3] two ways to compute R(x, v n ), 

reconstruction via primitive function (RP) and reconstruction via deconvolution (RD) were 
described. We only review RP here since a knowledge of it is enough for us to introduce our 
slope modification method. 

Notice that, for any fixed integer to, {p,}^, 0 = Vjh}?l iQ is a sequence of pointwise 

values of the primitive function p(x) = /* u(x)dx of u(x) at x, o _i ,x f - . i ,x, + s, • • • . A 
natural way of getting a polynomial which approximates u in the cell (xj_i, x J+ i ) is, there- 
fore, to interpolate p(x) at r + 1 consequent points x l(j) _i, x i( j )+ i, .. ., x,^ +r+ i, including 
Xj_ i and x j+ i. The derivatives of the interpolating polynomial then give the reconstruction 
and its derivatives. One has r degrees of freedom of choosing i(j), i.e., the stencil. It is the 
way of choosing the stencil that distinguishes the ENO schemes from others. 

For simplicity, we only give the hierarchical algorithm for evaluating i(j). The goal is to 
find the “smoothest” stencil which includes Xy _ ^ and x J+ ^. The smoothness of the stencil 
is somehow measured by the absolute value of the divided differences p[xj, x J+1 , • • • , x; +fc ] 
which are defined inductively by 

(2.2) p[xj] = p(ij) j = -oo, • • • , +oo, 

and 

p[xj,Xj+ i,--*,X j+jt ] = (p[x i+1 , • • • , X : +k\ - p[Xj, • • • , x J+fc _ 1 ]/(x j+fc - Xj), 

(2-3) k = 1,2, • ■ • ,r 

i = — oo, • • • , oo 

We describe the algorithm evaluating i(j) in the flavor of Fortran language as follows: 


5 


1) i(i) = j- 

2) For k = 1,2, ••• ,r, if 


then 




*0') = *(i) - 1- 

Havmg computed i(j), one then gets a polynomial pfix,v) by Newton’s form of interpo- 
lation. From this the reconstruction R(x, v) and its derivatives at x = xj are given by 
d‘ d l+1 

-j-jR(x,v)\ xmX) = -——-Pj(x,v) la— ^ , / = 0, 1,* **,r — 1. 


§3 The Slope Modification Method for Scalar Conservation Laws 


Denote by 6^ the jump of R(x, u) at the cell interface x y i, i.e., 


(3.1) 


~ +0, u) — R(xj_^ — 0,u). 


We introduce the following slope modification algorithm: 

Algorithm 3.1. The modified reconstruction R(x,v) in is a (r-l)-th order 

polynomial which is defined by 


(3.2) 


d l - d l 

^R(x,v) |*=x,= —R(x,v) |*= xy , 


l<r-l ,l± 1, 


( 3 - 3 ) ^R{x,v) 1 ^ 

where the slope modifier dSj is given by 


£*»■»> 


\x=x } +dSj/h, 


0.4) 


dS : = 2m(a J m(<J i _ i ,5 i+| ), m (t; i+1 - R] +} .,R +_ , - t^)). 


here {a^} is a sequence of positive numbers, the function m(x, y) is defined by 


(3.5) 


m(x,y) 


0, xy < 0, 

k m u i (|®|i |)sgn (x), otherwise, 
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and Rf+^ denotes /2(xj + ^ ± 0, v). 

To see the effect of the algorithm, let us apply it to the UNO scheme introduced in 
[1]. The UNO scheme is based on a non-oscillatory piecewise linear reconstruction. More 
precisely, if v = {ujjjl.,*, are the cell averages of a function it, the UNO reconstruction 
R(x, v ) can be written as 

( 3 - 6 ) R(X, V ) = Vj + Sj - h Xj , < x < x j+ ^. 

The reconstruction is non-oscillatory in the sense that the number of the extrema of R( ■ , v) 
is no more than that of v. We recall that R(-,v) satisfies(see [2]) 

(3-7) - Uj-i) > 0, 

(3.8) 6-_ ^ = 0, if Vj — vj-i = 0, 

(3.9) Sjxaivj+i - vj, Vj - vj _ t ) > 0. 

Recall also that (see [1]), if the numerical scheme for scalar conservation law (1.1) is 

(3.!0) »7 +1 = «; - A(/. +i - f^), 

a modified scheme adding artificial compression to (3.10) would be 

< 3 - n ) «r = »; - *(/** - /j-i). 

where the flux modifier 

( 3-12 ) 9j+$ = fj+% ~ /,+ ». 

obeys 

(3.13) <7j+^Au" > 0. 

Consider the UNO scheme with numerical flux(see [2]) 


(3.14)/,,, = | 

| f(v") + 0.5« i+i (1 - X^)Sf/[l + A(S /+ j - «,_*)] 

if S+3. > 0, 

where 

{ f( v j+ 1) ~ 0.5a j+ ^(l + Aa j+ |)S" +1 /[l + A(a j+ 3 - a j+ ^)] 

Ol 

+ 

A 

o 

(3.15) 

= (/(«* 1) - fWMvj+i ~ vj). 
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In the remaining part of this section, by ‘ the UNO scheme we always mean the scheme 
with this flux. 

Applying the slope modification method to the UNO scheme, one gets a scheme whose 
numerical flux has the same expression as (3.14) except that Sj is replaced by Sj for all j. 
From (3.4), (3.7), (3.8) and (3.13) one can see that under the CFL type conditions 

C 3 - 1 ^) |Aa i+ ^| < 1 

and 

( 3 ‘ 17 ) A(a i+ ^-5 i _^) > -1, 

algorithm 3.1. applied to the UNO scheme has the effect of adding artificial compression to 
the UNO scheme. 

Remark 3.2. For general high order ENO schemes, (3.7), (3.8) and (3.9) are not 
necessarily true. Nevertheless, the numerical experiments in [4] (see the figures of the recon- 
structions at the end of [4]) show that they are essentially true. Furthermore, the numerical 
results we are going to report in §5. show that the slope modifier with suitable chosen oij 

works well for all the cases that we have tested, although a rigorous analysis is not available 
now. 

Next, we discuss the stability of our method. Since the modified reconstruction no longer 
obeys (3.7), (3.8) and (3.9), the algorithm does introduce oscillations to the reconstruction. 
However, due to the cell average process, it seldom introduces oscillations to the solution at 
the next time level. In fact, we have 

Theorem 3.3. If \\a j+ ± - a^\ < 1 and A|5 j+ i | < 1 hold for all j, then the modified 
scheme obtained by applying algorithm 3.1. to the UNO scheme is non- oscillatory. 

In order to prove the theorem, we first describe the procedure for deriving the modified 
UNO scheme. This procedure is a direct copy of that for the UNO scheme. Having obtained 
the reconstruction R(x,v n ) and the corresponding modified reconstruction R(x,v n ), one 
derive pointwise approximate solutions w(x) and u>(ar) at time < n+1 as follows: both w(x) 
and ti>(x) are piecewise linear with possible discontinuities at 


(3.18) 


X }+\ ~ x i+ j + a j+± T 
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w ( X j+$ ± 0) = R(x j+ i ± 0, u"), 


such that 
(3.19) 

and 

(3-20) ut(X j+ i ±0) = R{x-^_ ± 0, v n ). 

We shall call these points of possible discontinuities d-points. The modified UNO scheme is 
then completed by defining 


(3.21) 


= U 




w(x)dx. 


For convenience, let us make some definitions and notations. 

Definition 3.4. A sequence of points x,-, Xj+i ,.. . ,Xj is called a strict maximum section 
of a cell average function v = and is denoted by S(v,i,j) if 


(3-22) u,_i < Vi = u i+1 = • ■ ■ = Vj > vj +l 

Denote by 5(u) the set of all such S(v ; i, j)'s. 

Definition 3.5. An interval [a, 6], a < b is called a strict maximum section of a 
piecewise linear function u and is denoted by s(u; a,b) if there exists a<5>0 such that 


u(x) < c = limsupu(£) 

holds for a — <5<x<aor6<x<6 + <5, and such that u(x) = c if a < x < b. Denote by 
s(u) the set of all such J(u;a,6)’s. 

Similarly, one can define strict minimum sections £L{v;i,j) and s(u;a,b) along with cor- 
responding sets £_(v) and s(u). 

Remark 3.6. In the Definition 3.5., if u is the UNO reconstruction of a cell average 
function v , and if a = 6, then one can make u(x) continuous at x = a by defining u(a) = c. 
The same observation holds for the strict minimum sections. 

Remark 3.7. If s(w; a, b) is a strict maximum section of w, then s(w; a, b) is also a strict 
maximum section of w. Furthermore, 

limsupu>(£) = limsupu;(£). 

a £-*a 

The same observation holds for the strict minimum sections. 
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We denote by N(u) the number of the strict maximum sections of a function u, by N(u) 
the number of its strict minimum sections and by N e (u) the number of all its strict extremum 
sections. Finally, we denote by u[a, 6] the average of a pointwise function u over (a, b). 

We have trivially, 

(3-23) N e (w) < N e (R(-,v n )). 

We will prove the theorem by showing that 

(3-24) iV e (u n+1 ) < N e (w). 

This inequality is implied by the following two inequalities 
(3-25) N(v n+1 ) < N(w) 

and 

(3.26) jV(u n+1 ) < W(to). 

Since the proofs of the two inequalities are same, it suffices to prove (3.25) which is equivalent 
to the following 

Claim. There is a 1-1 mapping 
(3-27) P:5(u n+1 ) s(u,). 

We need the following 

Lemma 3.8. Under the conditions of the theorem 3.3, If w increases(decreases) in 
(Xj_^,x i+ 3), then 

( 3 - 28 ) wgi 1 > (<)v” +1 . 

Proof. Assume that u;(x) monotonically increases. By (3.21), 

(3-29) vffi - u" +1 = i J_ h/2 (w(x j+1 + 0 - w( X] + t))d£. 

It suffices to show that 

(3-30) w(x j+ i + 0 > w(xj + £), for - h/2 < £ < h/2. 
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Notice that both Xj + £ and Xj± { + £ are in x J+ a). If ( Xj + £,xj+i + £) contains no 

d-point, inequality (3.30) is trivially true. If (xj + £, x j+1 + £) contains more than one d- 
point, say, X k+ ^, X k+ 3 , . . . , X k+m+ ^, then by algorithm 3.1. and the monotonicity of tu(x), 
w(xj + £) — u ?+i ^ + 0- Finally, suppose that ( xj + £, x J+1 4- £) contains only one 

d-point, say, X k+k . Under the conditions of theorem 3.3., the inequality 

0 < X i+\ ~ X i-\ < 2h 

holds for all j. Therefore, if we extend w(x) in X fc+ ^) linearly to -oo, extend w(x) 

( X k+ ^ , X k+ 3 ) linearly to +oo, and denote the extension by 1U ( x ) , then, 

(3.31) W(X k+i + 0) > W {XkH - A), 
and 

(3.32) W(X iH + A) > W(X kH - 0). 

The linearity of W{x) in each half-line devided by X k+ ^ implies that for the above <f, there 
is a 0, satisfying 0 < 0 < 1, such that 

(3.33) w( Xj + 0 = 0W(X k+k - h) + (1 - 9)W(X k+} . - 0). 
and 

(3.34) w(ij+i + 0 = 0W(X k+ ^ + 0) + (1 - 0)W(X k+ 1 + h). 

Clearly, (3.31)— (3.34) imply (3.30). The lemma is then proven. 

Proof of theorem 3.3. We define a mapping P from S'(u n+1 ) into s(u;) as follows. 

Given S(v n+1 ;i,j) £ S{v n+1 ) , one of the following five cases must occur. In each case, we 
will determine the image P(5(u n+1 ; i, j)) £ s(u>), and prove that the mapping is well-defined. 
Case 1. There is at least one strict maximum section s(w;a,b ) £ s(w ) such that 

[a, 6] n [x,_|,x i+ |] ^ <t>. 

and 

limsupin(x) > u” +1 . 
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In this case, we pick up this strict maximum section to be P(5(u n+1 ; i,j)). Obviously, 

(3.35) x,_ 3 < a < b < x J+ 3 

holds. 

. Case 2. The following two conditions hold: 

1) . There is no strict maximum section s(w;a,b ) £ s(tw) which satisfies the conditions of 
the case 1). 

2) . w(x) = u, n+1 , for x,_j. < x < x j+ ^ 

Let (c, d) be the largest interval containing (x,_^,x i+ |) such that w{x) is constant in it. 
Since (c, d ) is not a strict maximum section of w(x), w(x) must either strictly increase in 
some right neighborhood of x j+ ^ or strictly decrease in some left neighborhood of x^j.. If 
the former is the case, according to the definition of S(v n+1 ;i,j ), there must be a left most 
s(w; a, b) such that 

d < a < b < X : , a . 

j-r a 

We then define P(5(u n+1 ; i, ;')) to be this s(u;; a, b). Otherwise we can pick up the right most 
one as P(5(u n+1 ; i, j)). In either case, (3.35) holds. 

Case 3. u>(x) increases but is not constant in (x,_^, x j+ ^). We take P(5(u n+1 ; i,j)) to 
be the left most s(u;;a,6) € s(w) such that a > x j+ ^. The well-defined ness is justified by 
above lemma, and (3.35) holds. 

Case 4. w(x) decreases but is not constant in (x f _^,x i+ ^). We take P(S(v n+1 ; i, j)) to 
be the right most s(w;a,b ) <= s( w) such that b < x f _|. The well-definedness is justified by 
the same lemma, and (3.35) holds. 

Case 5. The following two conditions hold: 

1) . There is no strict maximum section s(tu;a,6) 6 s(w) which satisfies the conditions of 
the case 1. 

2) . w(s) is not monotone in (x,_^,x J+ |). 

In this case, there must be a £ in (x,_^,x j+ ^) which satisfies the following requirements: 
a), w = w in a neighborhood of f. b). u>(£) = w(£) < v? +1 and c). Either 

(3.36) u?(x,_^,(] < uf +1 
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or 


*(**-^€1 > u . n+1 


(3.37) 

is true. If (3.36) is the case, we define P(5(u n+1 ; i,j)) to be the left most s(w; a, b) 6 s(w) 
such that a > x j+ Otherwise, we define P(5(u n+1 ; i,j)) to be the right most s(w;a,b) € 
s(u>) such that b < x ; _|,. To see that such 5(it»; a, b) exists, let us assume that (3.36) holds, 
then we have 

(3.38) u>[£,x j+ ^] > t;" +1 . 

Suppose that k is the integer satisfying 

x k+\ ^ x j+\ > X k _i. 

If Sg < 0, then 

M = sup w(x) = sup w(x) > v" +1 . 

*6({,x j+ p x€(«.x j+ p 

Assume that w(x) has been made continuous at those d-points mentioned in the remark 
3.6.. Then, M will be attained by w(x ) in (£,x J+ ^). There is thus a strict maximum section 
s(w; a, 6) such that f < a < b < x . . i and 

limsupu>(x) > v? +1 . 

x—*a v y 1 

This contradicts the condition 1). Therefore one must have Sg > 0. 

Using a similar argument one can show that vg > vg +1 . Hence, if w(x) monotonically 
increases in (Xj+%, x j+ 3 ), then 

(3.39) wgi 1 >vg> v] +1 , 

which contradicts the definition of S(v n+1 ;i,j). This implies that there must be a s(w;a, b) 
such that 

x J+ 3 > b > a > x j+}r 

Similarly, if (3.37) holds, then, there must be a s(w;a,b) such that 

x,_| < a < b < x,_|. 
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Finally, we have to show that P is 1-1. Let $(v n+1 ; i x ,j x ), and S(v n+1 ;i 2 ,j 2 ) be two 
adjacent strict maximum sections of v n+1 , and 

P(S(v n+1 \i k ,j k )) = s(w;a k ,b k ),k = 1,2 
From the argument above, we see that if 

P(5(u n+1 ;z,;)) = J(u;; a, 6), 

then 

(3-40) x i-\ < a < b < Xj + 3 . 

Therefore, if a x = a 2 and b\ = b?, the only possibility would be 


and 


*2 — ji + 2 


x h+h < a l < ^1 < x h _i 


Now, one can see easily that either 


u " +1 > ?) n+1 
u ji+i - v h ’ 

which contradicts the definition of 5(u n+1 ; i\,j x ), or 


> < +1 

which contradicts the definition of 5(u n+1 ; i 2 ,j 2 ). The claim and, therefore, the theorem is 
proven. 

Next, we consider the order of accuracy of our method. For linear equation 

«< + u x = 0, 


we have 

Theorem 3.9. Suppose the aj’s in algorithm 3.1. are uniformly bounded, then at the 
worst, this method lowers the accuracy of the original scheme by one order. If in addition, 


a j + 1 ~ <*j = 0(h), 
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then the method applied to the ENO schemes keeps the order of the accuracy. For the choices 
of the parameters we are going to suggest, this last condition is valid if there are no zeros of 
the first r — 1 derivatives of the true solution u near Xj 

Proof. It suffices to recall that, at the cell interface x J+ |., the jump <5 J+ 1 . of the recon- 
struction satisfies 

(3.41) = c j+i h' + 0(h’+') 
and 

(3.42) Ac j+k = 0(h) 
if 

(jj) ,u ( z ) ^ °’ 1 ~ 1- 

We believe that the conclusion of above theorem is also true for the fully nonlinear 
problems, but do not have a proof yet. 

Next, we discuss the choice of the afis. Hereafter we shall call them the SM(slope 
modification) coefficients. Unfortunately, we have not been able to give a universal method 
for determining the afis. The recommendations we are going to give come from our numerical 
experiments. We found that, for linear problems, the smearing effect is essentially eliminated 
when the afis are greater or equal to 1.9 - 2.3. One usually gets satisfactory results by letting 
otj be equal to or slightly larger than the smallest number capable of eliminating the smearing 
effect. 


For the problems with rich smooth structures, the uniform choice of aj as suggested 
above could run into trouble. Although the essentially non-oscillatory property is generally 
kept, some artificial contact discontinuities may damage the smooth structure. 

One way of avoiding this is to use a discontinuity detector. A well known one is 

A 2 !!,'-! 


(3.43) 

Then aj can be written as 

(3.44) 


Mi = 


|Auj_i| + |A«j| 


aj = cpj 


In most computations we have performed, c = 33. 
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Remark 3.10. Usually, (3.43) and (3.44) give satisfactory results. However, for the local 
CFL number far away from ±0.5, the compression effect is not balanced between the head 
and the tail of a discontinuity. Further improvement can be made by multiplying the coeffi- 
cients (3.44) with a balance factor. For details of the balance factor, see §5. 

§4 Application to the ENO schemes for the Euler equations of gas dynamics 
In this section, we apply the slope modifying method to the ENO schemes for the Euler 
equations of gas dynamics for a polytropic gas. The application to the general systems of 
hyperbolic conservation laws follows immediately. 

For a polytropic gas, the governing equations are, as in [4], 


(4.1) 

u t + f{ U )x = 0 

(4.2) 

u = ( p,m,E) T 

(4.3) 

f(u) = q u + (0,P,qP) T 

(4.4) 

P = (7 -1)(£- 


Here p, q, P, and E are the density, velocity, pressure and total energy, respectively; m — pq 
is the momentum and 7 is the ratio of specific heats. 

The eigenvalues of the Jacobian matrix are 


( 4 - 5 ) a x (ti) = q - c, a 2 (u) = q, a 3 (u) = q + c 

where c = (~fP/ p)i is the sound speed. 

The corresponding right-eigenvectors are 




1 ] 


1 ] 


1 ) 

(4.6) 

r i( u ) = 

q-c 

{ H ~qc j 

,r 2 (u) = 


, r 3(«) = 

q + c 

\ H + qc t 


H=(E + P)/p = c>/( 7 -l) + i , 2 


here 

(4.7) 

is the enthalpy. 
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The corresponding left-eigenvectors bi-orthonormal to (4.6) are 



h(u) = 

j(f >2 + q/c, —biq - 

l/c,6i) 

(4.8) ■ 

h(u) = 

(1 - b 2 ,b 1 q,-b l ) 



k w = 

?(&2 ~ q/c,-b 1 q + 

1/cA) 

where 




(4.9) 


&i = (7 ~ 1 )/<? 


(4.10) 


h = 



To avoid too many collisions of the discontinuities which damage the advantage of the 
ENO reconstructions, ENO uses the characteristic reconstruction. Numerical experiments 
demonstrate that our slope modifier should also be applied in characteristic variables to get 
rid of some, although minor, spurious oscillations. In addition, since our goal is to sharpen 
the contact discontinuities which are the 2-waves, we only have to use it in the second field. 
We now introduce two algorithms to modify the r-th order characteristic reconstruction. 
Algorithm 4.1. For each j, 

1) . Compute the locally defined characteristic variables 

™H v j) = h( v")v? fori = j-r,...,j + r if k = 1,3, 

for i = j - r - l,...,j + r + 1 if k = 2. 

2) . Apply the scalar reconstruction algorithm to each of the locally defined characteristic 
variable in ( 4 . 11). The result is 


(4.11) 


(4.12) 


Ql 

bj.i ~ ™ ))!***>» / = 0, 1, . . . , r — 1; A; = 1, 3, 


& 


and 

( 4 - 13 ) ^i?(a;;u;V;))U=x J+m , / = 0,l,...,r — l;m= -1,0,1. 

When m = 0, / ^ 1, we denote (4-13) by b 2 t . 

3) . Apply our method, i.e., add the slope modifier to (13) to get b 2 - x . 

4) . Transform back to physical variables: 

( 4 - 14 ) -R(z; u n ) = £ b hl (x - Xj) l /l\ x-_ 1 < X < x j+ i 

i * 2 


/=0 
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where 

(4.15) 


3 


hi = b j,l r k- 


k=l 


Algorithm 4.2. 1). For each j, do the ENO characteristic reconstruction as usual to 
getR(x,v*) for Xj _ h <x<x J+k . Store l 2 (v*), r 2 (vf), uj(vf), w j +1 (v]), and 

b h- 


2) . For each j , compute = R( Xj+ ^ ± 0,u"), = R+ +1 _ - i?" 

3) . For each j, 

i). compute 


R + (w) = kiv^Rj^, 
R~{w) = h(vj)R ^_ , , 

j y 

= h{vi)6 i+k . 


8 (w) = l 2 {vj)6-_^\ 

ii). using wj.tivf), wj(v]), tifi j+l (vj), R+(w), R~(w ), 8 + (w ), and 8~(w) to get 
the slope modifier db 2 - lf 

Hi), compute bji = hi + dblMvf). 


§5 Numerical results 

Since our purpose in designing the slope modification method is to improve the per- 
formance of the ENO schemes for the solutions which contain contact discontinuities, the 
numerical experiments we will present in this section are either on linear equations with dis- 
continuous initial values or on the Euler equations for gas dynamics whose solutions contain 
contact discontinuities. All the examples are standard problems widely used in the literature 
to test various schemes for the hyperbolic conservation laws. 

With these examples, in addition to showing the performance of our method, we also try 
to give some suggestions on how to choose the parameters in our algorithms since we have 
not been able to establish any theoretically meaningful rule to determine these parameters. 
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Example 5.1. We apply our method to ID linear model problem 


(5.1) u t + u x = 0, 1 < x < 1 

(5.2) u(x,0) = u o (x) 
with periodic boundary condition. 

We test our method for this problem with each of the following four initial conditions. 
Divide the interval [—1,1] into 100 cells of equal size. The center of the j-th cell is Xj = 
(j — 50 .5) h, where h = 1/50. The four initial conditions are 

1, 35 < j < 65, 

0, otherwise, 

50)/15] 2 ) 1/2 , 35 < j < 65, 

otherwise, 


(5.5) 


v 0 _ e — 300(x ; -0.5) 2 


and 






—yj sin( f^ry 2 ), 

-1 < 

(5.6) 

- < 

v j+ 25 (mod 100) — < 

| sin(2xt/j)|, 

l»il < 



2 Vj - 1 - sin(3iry i )/6, 

5 < Vj < 1 


In the last condition, t/j = x ; — 0.01. Essentially, the first 3 conditions are those used by 
Zalesak [15] (see also the references therein) and the last one is that used by Harten et al 
[4]and [9] 

The numerical solutions for the initial condition (5.3) are displayed in the Figures 1 - 
19. In the computations for the Figures 1-5, the 2nd order ENO scheme is used. The CFL 
number is fixed at 0.8 and the SM coefficients a : ’s in (3.4) are chosen to be independent 
of j and n. We increase ctj = a from 0 to 10 and, for each choice of a, run the program 
twice for 250 timesteps and 1250 timesteps respectively. The figures show clearly that, for 
small a, the discontinuities are smeared more and more as the the number of the time steps 


(5.3) 


”?=< 


(5.4) 


”?=< 


(1 - [U - 

0 , 
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increases. However when the coefficient reach about 1. 9-2.3, the smearing appears to be 
essentially eliminated. This is indicated by the coincidence of the graphs of the numerical 
solutions obtained by running the program over the two different numbers of time steps. For 
a given ENO scheme and a fixed CFL number, we call the smallest a which is capable of 
eliminating the smearing effect the critical value of the SM coefficient with respect to the 
scheme and the CFL number. 

Remark 5.1. A quantity equal to or slightly larger than the superemum of the critical 
value over the range [0,1] of the CFL number is a good candidate for the choice of the 
MS coefficient. For the solutions which lack smooth structure we can choose a larger SM 

coefficient. Otherwise, we should choose one near or equal to the superemum of the critical 
value. 


Remark 5.2. For the solution dependent SM coefficient (3.44), one can similarly “define” 

the critical value for the parameter c. The recommondations in the last remark also applies 
here. 

The Figures 6-13 show the effect of CFL numbers on the critical values. We do the same 
computations as for the Figure 1-5 with different CFL numbers. The SM coefficients are 
Q = 1*9 3-nd a. = 2.3. We found that the effect is minor. 


Notice that although the CFL numbers have little effect on the critical value of the SM 
coefficients, it does play a role in the profiles of the solutions over the linear discontinuities. 
The sharpening effect is not well balanced, i.e., the steepness at the heads of the disconti- 
nuities is different from that at the tails when the CFL numbers are far away from ±0.5. 


Looking at the Figure 1 carefully, one sees that this property of nonsymmetry already exists 
with the original ENO schemes. 

For solutions with rich smooth structures and/or for multi-dimensional problems, when 
our slope modification method is applied and the resulted programs are run for huge number 
of timesteps, this phenomenon becomes more serious. To overcome this trouble, observe that 
at the head of a discontinuity, is very large, while at the tail, it is very small. We 

hence define the following balance operator 

I Aj— il>|6(A J -0.5Sgn(A J )) 


(5.7) 


Wi,6) = 


AjV 
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where 6 is a positive parameter to be chosen, \j = a.^ is the local CFL number, and a : is 
the characteristic speed. We will see the effect of this operator in the examples displayed in 
the Figures 26 - 29. 

The Figures 12-19 show the effect of the order of accuracy on the the critical values of 
the SM coefficients. In these figures, r-0 denotes the “ r-th order ” ENO schemes, while 
r-1 denotes the r-th order ENO schemes, both in the sense of [4] . We find that the effect 
is also minor. This is interesting since the jumps of a higher order reconstruction at the cell 
interfaces over the transitions of the discontinuities are much smaller than those of a lower 
order reconstruction (see, again, the figures at the end of [4]). 

The fact that the critical value of the SM coefficients is relatively independent of rather 
than inversely proportional to the jumps for different orders of the schemes is consistent to 
the fact that the higher order ENO schemes smear the contact discontinuities less then the 
lower order ones do. 

The results for the initial condition (5.4) are displayed in the Figures 20 - 29. In the 
Figures 20 - 25, we use the SM coefficient 1.9. The improvement is apparent. The Figures 
26 — 29 test the effect of the balance operator (5.7). This time we run our program of the 
2nd order ENO for 6000 timesteps. For the Figures 26 and 27, we use the constant SM 
coefficients 1.9 while for the Figures 28 and 29, we use (3.44) with c = 33. The results 
obtained by using the SM coefficients multiplied by the balance operator with b = 4.3 are 
displayed in the Figures 26 and 28. The Figures 27 and 29 show the corresponding results 
without applying the operator. 

The Figures 30-31 displays the numerical results for initial condition (5.5); the Figures 
32-33 are the results for initial condition (5.6). Again one can see apparent improvements 
of our method in the treatment of the cusps, jumps or both. 

Example 5.2 We apply our method to the Riemann problems for the Euler quations of 
gas dynamics (4.1) with following two sets of initial conditions known as the Sod problem 
and the Lax problem respectively: 

^ {PLi 9l? F l ) (1? 9» 1)> 

(PR,qR,PR) = (0.125,0,0.10) 
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and 


( 5 9 J (pL, qu Pl) = (0.445, 0.698, 3.528); 

( Pr , Pr) = (0.5, 0, 0.571) 

Both problems are solved by the “4-th order” ENO schemes with and without the slope 
modification. In all the computations, we use 100 equally sized space cells and the CFL 
number 0.8. The Figures 34 and 35 are the numerical results for Sod problem after 60 time 
steps; the Figures 36 and 37 are those for Lax problem after 85 timesteps. Here we apply 
the alghorithm 4.1 with the SM coefficients (3.44) such that c = 33. Notice that with the 
slope modification method the computed contact discontinuities are of the same quality as 
the discontinuities captured by the same method in the linear model problems above. 

Example 5.3 (The blast waves problem). Here we consider the equations (4.1) with the 
following initial condition 


(5.10) 


u l, 


u(x, 0) — < um, 


{ 


0 < x < 0.1 
0.1 < x < 0.9 
0.9 < x < 1 


where 

( 511 ) PL = Pm = Pr = 1 , qL = qM = qR = 0 , 

^ = 10 3 ,Pa/ = 10- 2 ,P* = 10 2 , 

and the two boundaries are assumed to be solid walls. See [14] for the details of the solution 
and the comparison of the performance of various schemes for this problem. Notice that 
the contact discontinuity results from the collision of the two strong shocks. It is extremely 
difficult to be “captured” by a shock capturing scheme in the Eulerian framework. In fact, no 
scheme tested in [14] captured the contact discontinuity well. To the auther’s knowledge, in 
the literature of the modern shock capturing technique, the only successful computations of 
this discontinuity so far have been performed by Harten’s subcell resolution(see [9] and [7]). 
In the computations for this problem, we use 200 cells and the CFL number 0.8. The results 
at t = 0.038 are displayed in the Figures 38, 39 and 40. The solid lines are the numerical 
solutions obtained by the 2nd order ENO scheme and the slope modification method with 
800 cells. The Figure 38 is the result of the 2nd order ENO scheme with the algorithm 4.2 
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and the Figure 39 is the result of the 4-th order ENO schemes with the algorithm 4.1 . For 
comparison, in the Figure 40 we display the numerical result of the 4-th order ENO without 
the slope modification. Notice that all the contact discontinuities are well captured by the 
ENO schemes with our method. Notice also that the two algorithm perform equally well. 

Example 5.4 To see the performance of our method for the problems that have some 
structure, we apply it to the Euler equations (4.1) with the following initial condition: 


(5.12) 


(P, <1,P) 


(3.857143, 2.629369, 10.33333), x < -4 
(1 + 0.1 sin5x,0, 1), x > — 4 


This is a model problem for shock/ turbulence interaction. See [12] for a linearized analysis 
of this problem, and [7] for a numerical result. We apply our slope modification method 
with the balance operator on the “4— th order” ENO schemes. The results are demonstrated 
in the Figures 41-34. The computations are performed with the CFL number 0.8. The solid 
line is the result with 800 cells. Comparing the result with that in [7] we can regard it as the 
exact solution. The circles in the Figures 41 and 42 are the results computed with 200 cells 
and 400 cells respectively. For comparison, the Figures 43 and 44 show the results computed 
with 200 cells and 400 cells respectively without the slope modification. 

Example 5.5 This is a preliminary result for 2D computations. Consider the 2D model 
problem 

(5-13) u f + u r + u v = 0, -1 < x, y < 1 

with the initial condition 


(5.14) 


«o(x,y) = < 


l*-y| < j 


[ 0, otherwise. 

In [10], Harten used this problem to test the performance of the ENO schemes with the 2D 
reconstructions via deconvolution. Now, we use the reconstructions via primitive functions 
dimension by dimension and, at the same time, apply the scalar slope modification algorithm 
with the balance operator. The CFL number is 0.8 x 0.8. The computations are performed 
with 20 x 20 cells. The results with the slope modification are displayed in the Figures 45-47, 
while those without in the Figures 48-50. The apparent improvements shown in these re- 
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suits and the results in [7] indicate that the pfesent method is promising in 2D computations. 
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Figure 44.1: shock/turbulence, 4-th order ENO without SM 
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Figure 44.3: shock/turbulence, 4-th order ENO without SM 
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Figure 45: 2nd order ENO, with artificial compression 



Figure 46: 3rd order ENO, with artificial compression 


Figure 47: 4th order ENO, with artificial compression 



Figure 48: 2nd order ENO, without artificial compression 


Figure 49: 3rd order ENO, without artificial compression 



Figure 50: 4th order ENO, without artificial compression 



